Objectifs de la formation

  1. Comprendre le séquençage Illumina et les fichiers FASTQ
  2. Manipuler des données de séquençage avec R
  3. Découvrir BARQUE, un pipeline de métabarcoding québécois
  4. Appliquer les concepts sur des exemples pratiques

Avant de commencer

Installation des librairies R requises

  1. Ouvrir RStudio
  2. Envoyer les commandes suivantes dans la console
# Installation des packages Bioconductor
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install(c("ShortRead", "Biostrings", "Rqc"))

Introduction

Comment le séquencage est effectué?

Crédit: J. Couillard

DOI: 10.1016/B978-0-12-802234-4.00002-1

Comment le séquencage est effectué?

Note

  • Chaque nucléotide à un fluorophore spécifique (4 canaux).
  • Certains appareils (ex. NextSeq) vont utiliser deux canaux car c’est plus rapide et moins coûteux.

DOI: 10.1016/B978-0-12-802234-4.00002-1

Le signal renvoyé par le laser est traduit en chromatographe pour chaque amplicon séquencé.

Comment le séquencage est effectué?

Exemple de chromatographe

Source: LabXchange.org

Comment le séquencage est effectué?

Exemple de chromatographe

Important

Le score est basé sur:

  • La séparation entre les clusters (pureté du signal)
  • Le rapport signal/bruit
  • La distribution de l’intensités des 4 canaux
  • La position dans le read (la qualité diminue généralement vers la fin)

Source: wikipedia.org

Fichiers de sortie du séquenceur Illumina

Exemple de nomenclature de fichier de sortie du séquenceur Illumina: sample1_1_L001_R1_001.fastq.gz

Signification Exemple
sample1 Identifiant de l’échantillon sample1, sample2
_1 Numéro de réplicat 1, 2, 3…
L001 Numéro de ligne / lane: L001-L008
R1/R2 Direction de lecture R1=avant, R2=arrière
001 Segment de fichier 001, 002…

Note

  • Séquençage paired-end = fichiers R1 + R2 pour chaque échantillon avec l’ensemble des amplicons.
  • Numéro de Ligne ou lane = Une flowcell Illumina est divisée en plusieurs voies physiques parallèles où le séquençage se déroule simultanément.

Structure du fichier FASTQ

Contient plusieurs lignes dont un ensemble de 4 lignes correspond à un amplicon.

system.file(package="ShortRead", "extdata/E-MTAB-1147/ERR127302_1_subset.fastq.gz") |> 
  gzfile() |> 
  readLines(n = 12)
 [1] "@ERR127302.8493430 HWI-EAS350_0441:1:34:16191:2123#0/1"                  
 [2] "GTCTGCTGTATCTGTGTCGGCTGTCTCGCGGGACATGAAGTCAATGAAGGCCTGGAATGTCACTACCCCCAG"
 [3] "+"                                                                       
 [4] "HHHHHHHHHHHHHHHHHHHHEBDBB?B:BBGG<DDAA?AABFEFBDBD@DDECEE3>:?;@@@>?=BAB?##"
 [5] "@ERR127302.21406531 HWI-EAS350_0441:1:88:9330:2587#0/1"                  
 [6] "CTAGGGCAATCTTTGCAGCAATGAATGCCAATGGGTAGCCAGTGGCTTTTGAGGCCAGAGCAGACCTTCGGG"
 [7] "+"                                                                       
 [8] "IIIIHIIIGIIIIIIIHIIIIEGBGHIIIIHGIIHIIIIIIIHIIIHIIIIIGIIIEGIIGBGE@DDGGGIG"
 [9] "@ERR127302.22173106 HWI-EAS350_0441:1:91:10434:14757#0/1"                
[10] "TGGGCTGTTCCTTCTCACTGTGGCCTGACTAAAACCCAGTGGCATTAAGAAAGAGTCACGTTTCCCAAGTCT"
[11] "+"                                                                       
[12] "GGHBHGBGGGHHHHDHHHHHHHHHFGHHHHHHHHHHHHHHHHHHHHHGHFHHHHHHHHHHHHHH8AGDGGG>"

Structure du fichier FASTQ

system.file(package="ShortRead", "extdata/E-MTAB-1147/ERR127302_1_subset.fastq.gz") |> 
  gzfile() |> 
  readLines(n = 4)
[1] "@ERR127302.8493430 HWI-EAS350_0441:1:34:16191:2123#0/1"                  
[2] "GTCTGCTGTATCTGTGTCGGCTGTCTCGCGGGACATGAAGTCAATGAAGGCCTGGAATGTCACTACCCCCAG"
[3] "+"                                                                       
[4] "HHHHHHHHHHHHHHHHHHHHEBDBB?B:BBGG<DDAA?AABFEFBDBD@DDECEE3>:?;@@@>?=BAB?##"
  1. Ligne [1] : Identifiant de séquence (commence par @)
  2. Ligne [2] : Séquence d’ADN (A, T, G, C) de l’amplicon
  3. Ligne [3] : Séparateur (+)
  4. Ligne [4] : Scores de qualité (un par base)

Ligne d’identification

Décomposition de l’identifiant

@SH00321:6:BWR98207-2813:1:1101:1065:1015 1:N:0:AACCATAGAA+GGCGAGATGG
  • SH00321 : ID de l’instrument
  • 6 : Numéro d’exécution
  • BWR98207-2813 : ID de la flowcell
  • 1:1101:1065:1015 : Tuile et coordonnées X:Y
  • 1:N:0 : Numéro de lecture, flag de filtre, bits de contrôle
  • AACCATAGAA+GGCGAGATGG : Codes-barres d’échantillon (double indexation)

Interprétation des scores de qualité

Exemple de chaine de caractères donnant le score de qualité (1 caractère par nucléotide): GGGGGGGGGGGGGGGGGGGGGGGG9GG-G9G9G9GGGG

Caractère ASCII Score Phred Taux d’erreur Précision
G 38 0,016% 99,984%
9 24 0,40% 99,60%
- 12 6,31% 93,69%

Tip

Règle générale : Q30 ou plus = haute qualité (99,9% de précision)

Scores de qualité du séquencage

Les scores de qualité utilisent l’encodage ASCII.

 Quality encoding: !"#$%&'()*+,-./0123456789:;<=>?@ABCDEFGHI
                   |         |         |         |         |
    Quality score: 0........10........20........30........40

Scores de qualité du séquencage

Avec R, comment convertir l’encodage ASCII en score de Phred et en probabilité d’obtenir une erreur (précision)?

# Convertir la chaîne de qualité en scores Phred
quality_string <- "GGGGGGG99GG-GGG"
quality_scores <- as.integer(charToRaw(quality_string)) - 33

# Afficher la correspondance
data.frame(
  Caractere = strsplit(quality_string, "")[[1]],
  Score_Phred = quality_scores,
  Precision = paste0(round((10^(-quality_scores/10)) * 100, 2), "%")
)

Scores de qualité du séquencage

   Caractere Score_Phred Precision
1          G          38     0.02%
2          G          38     0.02%
3          G          38     0.02%
4          G          38     0.02%
5          G          38     0.02%
6          G          38     0.02%
7          G          38     0.02%
8          9          24      0.4%
9          9          24      0.4%
10         G          38     0.02%
11         G          38     0.02%
12         -          12     6.31%
13         G          38     0.02%
14         G          38     0.02%
15         G          38     0.02%

Lire un fichier FASTQ avec R

# Lire le fichier FASTQ
fq <- system.file(package="ShortRead", "extdata/E-MTAB-1147/ERR127302_1_subset.fastq.gz") |> ShortRead::readFastq()

# Examiner la structure
fq
class: ShortReadQ
length: 20000 reads; width: 72 cycles

Extraire les informations principales

# Nombre d'amplicons
length(fq)
[1] 20000
# Extraire les amplicons
sequences <- ShortRead::sread(fq)

# Extraire les identifiants
ids <- ShortRead::id(fq)

# Extraire les scores de qualité
qualities <- Biostrings::quality(fq)

# Premier identifiant et séquence
ids[1]
BStringSet object of length 1:
    width seq
[1]    53 ERR127302.8493430 HWI-EAS350_0441:1:34:16191:2123#0/1
sequences[1]
DNAStringSet object of length 1:
    width seq
[1]    72 GTCTGCTGTATCTGTGTCGGCTGTCTCGCGGGAC...GTCAATGAAGGCCTGGAATGTCACTACCCCCAG

Résumé statistique d’un fichier FASTQ

# Longueur des séquences
seq_lengths <- Biostrings::width(sequences)
summary(seq_lengths)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
     72      72      72      72      72      72 
# Distribution des scores de qualité
qual_matrix <- as(qualities, "matrix")
mean_quality <- rowMeans(qual_matrix)
summary(mean_quality)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
   2.00   33.53   37.35   34.84   38.85   39.99 
# Composition en bases
Biostrings::alphabetFrequency(sequences, as.prob = TRUE)[1:4, 1:4]
             A         C         G         T
[1,] 0.1944444 0.2638889 0.2916667 0.2500000
[2,] 0.2361111 0.2222222 0.3194444 0.2222222
[3,] 0.2361111 0.2638889 0.2222222 0.2777778
[4,] 0.2500000 0.2638889 0.1666667 0.3194444

Assurance qualité avec la librairie Rqc

fq_files <- system.file(package="ShortRead", "extdata/E-MTAB-1147/") |> list.files(full.names = TRUE)
qa <- Rqc::rqcQA(fq_files)
Rqc::rqcCycleQualityPlot(qa)

Analyser la qualité du séquencage

Pourquoi la qualité des lectures décroit au fur et à mesure des itérations?

  1. Dégradation des fluorophores
  2. Accumulation d’erreurs de phase, soit la désynchronisation des clusters
  3. Épuisement des réactifs, compétition pour les nucléotides par exemple
  4. Dommages à l’ADN par stress chimique par exemple
  5. Accumulation de sous-produits (ex. Résidus de réaction s’accumulent)

Pratique - Analysez vos propres données

Note

Exercice simple: Utilisez le code R que nous venons de voir pour analyser vos propres fichiers FASTQ

  1. Créer un dossier d’analyse
  2. Placez vos fichiers FASTQ dans sous-dossier data/
  3. Exécutez les chunks (le code des diapos) pour analyser la qualité
  4. Comparez les résultats avec les exemples montrés
  5. Bonus: Explorer les autres graphiques dans le package Rqc

Tip

Si vous n’avez pas de données, vous pouvez télécharger l’exemple ici. Ce sont 10 échantillons de métabarcoding mitofish-12S

10:00

Points clés à retenir

  • FASTQ = format standard pour les données de séquençage
  • 4 lignes par lecture d’amplicon : ID, lecture de l’amplicon, séparateur, qualité de la lecture
  • Paired-end : les fichiers R1 et R2 contiennent les lecture d’amplicons à apparier
  • Scores de qualité : Plus élevé = meilleur (viser Q30+)
  • La qualité se dégrade vers la fin des lectures (normal)

Qu’est ce que BARQUE?

Un large éventail d’outils en bioinformatique

Hakimzadeh et al, 2024

Qu’est ce que BARQUE?

  • Programme écrit par Éric Normandeau (IBIS, ULaval)
  • Disponible sur Github depuis 2017
  • C’est un pipeline de traitement de données de séquences spécialisé sur le metabarcoding
  • Il agit comme chef d’orchestre d’un ensemble de programmes

Qu’est ce que BARQUE?

Creative Commons
  • FLASH v1.2.11+: fusion des lectures paired-end
  • VSEARCH v2.14.2+: outil d’analyse de séquences (version obligatoire)
  • Trimmomatic: outil de filtrage et découpage des lectures
  • NCBI BLAST: annotation des séquences non-correspondantes
  • Un ensemble de programmes UNIX comme cut, sort, uniq etc. qui manipulent des chaînes de caractères ou des nombres
  • R & Python: Visualisation et synthèse des résultats

Les grandes étapes de BARQUE

Vue d’ensemble du pipeline de métabarcoding

Étape 1: Validation du projet

  • Vérification des fichiers d’entrée
  • Validation de la nomenclature : SampleID_*_R1_001.fastq.gz
  • Vérification des logiciels
  • Présence de la base de données correspondante à l’amorce

Étape 2: Nettoyage des données

  • 2a. Nettoyage des lectures avec Trimmomatic
    • Suppression des bases de faible qualité
    • Découpage selon les paramètres définis
  • 2b. Retrait des amplicons non-appariés

Étape 3: Appariement des lectures

  • Appariement des amplicons (R1 + R2) avec FLASH
  • Création d’une lecture consensus pour chaque amplicon
  • Amélioration de la qualité globale

Étape 4: Séparation des amplicons

Pour chaque lecture fusionnée:

  • Recherche du primer forward au début
  • Recherche du primer reverse à la fin (complément inverse)
  • Calcul de la distance avec les primers de référence
  • Vérification de la taille d’amplicon attendue
  • Sortie: Un fichier par type d’amplicon

Étape 5: Retrait des chimères

  • Détection des séquences chimériques avec VSEARCH
  • Suppression des artefacts de PCR
  • Conservation des séquences authentiques

Étape 6: Assignement taxonomique

  • Attribution des espèces via bases de données correspondantes à l’amorce
  • Utilisation de BOLD (COI), MitoFish (12S), etc.
  • Génération des tables d’occurrence par espèce

Télécharger BARQUE

Pratique

Télécharger BARQUE en cliquant sur le lien

Ou en rendant à l’adresse suivante: https://github.com/enormandeau/barque

Les étapes en tant qu’utilisateur

  1. On place les fichiers fastq.gz avec la bonne nomenclature dans le dossier 04_data
  2. On configure les amorces (P3 et P5) utilisées avec le fichier 02_info/primers.csv
  3. On ajuste les paramètres du pipeline avec le fichier 02_info/barque_config.sh
  4. On s’assure que la base de données BOLD est téléchargé si on travaille avec les amorces COI.
  5. On éxecute le pipeline

Tip

Si vous voulez réutiliser la configuration du pipeline pour d’autres fois, c’est une bonne pratique de dupliquer le fichier et de le sauvegarder sous un autre nom.

Étape 1 - Placer les fichiers dans 04_data

On peut faire tout de suite l’étape 1 avec vos données.

Tip

Si vous avez pas de données, vous pouvez télécharger celles fournies par Étienne Normandeau, en cliquant ici. Ce sont 10 échantillons de métabarcoding mitofish-12S, chacun.

Étape 2: Configuration des amorces

On configure les amorces (P3 et P5) utilisées avec le fichier 02_info/primers.csv

TODO: Montre le fichier

Étape 3: Ajustements des paramètres du pipeline

On ajuste les paramètres du pipeline avec le fichier 02_info/barque_config.sh

TODO: Montre le fichier

Espèces intégré à la base de données par type d’amorces

TODO: Passer en revue les bases de données déjà dans BARQUE

COI

TODO: Aller chercher la base de données NCBI pour les primers COI

Pourquoi utiliser docker

BARQUE requière plusieurs programmes qui ne sont pas tous compatibles avec l’environnement Windows

Hakimzadeh et al, 2024

Comment docker fonctionne?

@Soroosh Nazem

Étape 1. Télécharger l’image docker

Étape 2. Executer l’image

Étape 3. Executer BARQUE dans l’image docker

Application Shiny

Exécuter BARQUE de manière embarqué

Application Shiny

Exécuter BARQUE de manière embarqué

Rapport: Obtenir le rapport de résultats